Global and local condensate and superfluid fractions of a few hard core Bosons in a 

combined harmonic optical cubic lattice 
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We explore the global and local condensate and superfluid (SF) fractions in a system of a few 
hard core (HC) bosons (TV — 8 and TV = 40) trapped inside a combined harmonic optical cubic 
lattice (CHOCL) at T — K. The condensate fraction (CF) is computed for individual lattice wells 
by separating the one-body density matrix (OBDM) of the whole system into components at the 
various lattice sites. Then each "lattice-site" component is diagonalized to hnd its eigenvalues. The 
eigenvalues are obtained by a method presented earlier [Dubois and Glyde, Phys. Rev. A 63, 023602 
(2001)]. The effects of interference between the condensates in the lattice wells on the CF in one well 
is also investigated. The SF fraction (SFF) is calculated for TV = 40 by using the diffusion formula of 
Pollock and Ceperley [Pollock and Ceperley, Phys. Rev. B 36, 8343 (1987)]. Our chief result is an 
opposing behavior of the global CF and SFF with increasing lattice wave vector k. In addition, the 
CF in a lattice well is enhanced by the interference with its neighbor wells beyond the result when 
the interference is neglected. The global SF is depleted with a rise of the repulsion between the 
bosons, yet at very strong interaction superfluidity is still present. The global CF remains almost 
constant with increasing HC repulsion. A reduction in the lattice dimension, i.e. an increase in the 
lattice wave vector, increases the local CF in each lattice well, but reduces the corresponding local 
SFF. At large HC repulsion, a coexisting SF-(vacuum)MI phase is established. 

PACS numbers: 67.85.-d,67.85.Hj,03.75.-b,03.75.Lm 
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I. INTRODUCTION 

The Bose-Einstein condensation (BEC) of bosons in 
optical lattices (OLs) has recently become a topic of great 
interest which motivated substantial work (lHl^. Inter- 
esting investigations included an experimental realization 
of the latter |l| , BECs in tight binding bands of OLs with 
different geometries @ , the Bose- Hubbard model (BHM) 
[i"2T - [l~5j , effects of the lattice dimension [ll| , and instabil- 
ities of BECs in moving two-dimensional (2D) OLs [H . 

A topic of importance is the measurement of the con- 
densate fraction (CF) in individual OL wells, for which 
only a few investigations have been reported 0, [13, EH ■ 
Investigations concentrated mostly on measuring the to- 
tal CF of the whole lattice system. For example, in a 2D 
lattice boson system, Spielman et al. [19( experimentally 
measured the CF of Rb atoms as a function of the lattice 
depth. Furthermore, Fang et al. [1C§ measured the CF of 
a 87 Rb gas released from an OL in a time-of-flight (TOF) 
experiment. 

Importantly, Chen and Wu [|| noted that most of the 
theoretical work focused only on BECs in ID OLs, and 
that the literature on 2D and three-dimensional (3D) OLs 
is still scarce. We were thus motivated to explore BECs 
in 3D OLs; the work of Brouzos et al. [H| provides also 
the chief motivation. 

According to Chen and Wu [4], a BEC confined in a 
strong ID lattice can be regarded as a chain of weakly 
coupled condensate islands trapped in the lattice wells. 
This idea has also been propagated by Shams and Glyde 
20] and there would be little tunneling between these 
islands, as Chen and Wu stated. Our upcoming formula- 



tion of the present problem is based on the latter thought. 
That is, we study the role of the interference between the 
condensate in one lattice well and the condensates in all- 
neighbor wells. In this regard, we were also motivated by 
the work of Baillie and Blackie [2l[ . 

Further, superfluidity in OLs has also been given con- 
siderable interest HEHH-ii]. Particularly the super- 
fluid (SF) to Mott-insulator (MI) transition [H, l2rj434| 
has been explored intensively ever since its first experi- 
mental realization (22j, as well as coexisting SF and MI 
domains of harmonically t rapped hard core (HC) and 
lattice boso ns flil . HH l33l . l35l. |36||. For example, Roth 
and Burnet [24| . used a twisted-boundary condition ap- 
proach to compute the superfluid fraction (SFF) in a one- 
dimensional (ID) lattice. Similarly, Hen and Rigol [26| 
applied a twist in the boundary conditions to evaluate 
the SFF of HC bosons in a superlattice. The twist in the 
boundary conditions can, however, be only applied to a 
homogeneous OL system. In addition, using the BHM, 
the effects of an external harmonic trap on the state dia- 
gram of lattice bosons have been explored by Rigol et al. 
14], specifically for the coexistence of MI and SF domains 
inside the system. The vanishing of the CF at large val- 
ues of the lattice depth indicated also the crossing from 
the SF to the MI domain 

The role of the HC diameter in the SF and conden- 
sate depletion of hard sphere (HS) bosons in OLs has, 
to the best of our knowledge, rarely been outlined be- 
fore, particularly in the SF-MI transition of a few-body 
system. So far, the investigations of the CF and SFF 
have been largely as functions of the lattice barrier (c.f. 
Shams and Glyde [2(| and Spielman et al. fl9|): in con- 
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trast we explore them as functions of the boson HC in- 
teractions. One of the few investigations on the role of 
the interactions in the condensate properties inside an 
OL was undertaken by Snoek et al. [371 ] and Ramanan 
et al. [33] for a ID Bose gas in an inhomogeneous OL. 
Snoek et al. explored the effect of interactions on the 
condensate properties of a Bose-Fermi mixture trapped 
in a 3D harmonic OL. Among their findings, was that 
the condensate is depleted with a rise of the Bose-Fermi 
repulsive interactions which are similar to the boson HC 
interactions. Whereas Snoek et al. conducted their cal- 
culations for the whole of their system, we explore the 
condensate at individual lattice sites and then sum all 
individual contributions. 

One important investigation most relevant to ours was 
presented by Brouzos et al. jl8j, who conducted studies 
on homogeneous few-particle bosonic systems in a ID 
multiwell trap. Some of our findings for a 3D OL are 
similar to theirs. Further, the latter authors articulated 
that exact studies of trapped bosonic systems "are par- 
ticular for a few number of particles" . Few particles have 
also been used, e.g. by Fang et al. [13] who did calcula- 
tions on a ID Bose gas in an OL with only ten particles 
in ten lattice sites. We shall devote a special section for 
connecting our results with the findings of Ref. [l8| later 
on below. 

In our work here, the CF and SFF of a few HC bosons 
trapped inside a combined harmonic optical cubic lat- 
tice (CHOCL) is studied at T = K. The cubic OL 
has 3x3x3 lattice sites, or lattice wells. Most impor- 
tantly, the nonspherical symmetry of the CHOCL system 
is dealt with which poses a difficulty for the computation 
of its global CF. This fact has been outlined recently 
by Astrakharchik and Krutitsky (38[, who presented a 
method for the calculation of the total CF of an inho- 
mogeneous Bose gas confined by an OL. Yet this work 
here proposes to evaluate the CF in each individual lat- 
tice well. Hence, another goal of our paper is to present 
methods for computing the CF and SFF in each well of 
the CHOCL system. The total CF and SFF are also 
computed. An important point to emphasize, is that we 
distinguish clearly between global and local CF and SFF. 
Motivated by the work of Xue et al. [l6| , we also decided 
to check the role of the lattice dimension on the prop- 
erties of the CF and SFF. In another attempt, we seek 
a MI state in our systems by going to large HC repul- 
sions and OL depths. We also make a clear distinction 
between the effects of repulsive forces with zero range, 
as encountered in the BHM, and repulsive forces with a 
nonzero range, as described by the HS Jastrow function, 
Eq.© below. A pair of bosons interacting solely by a 
repulsive delta function can still be brought close to each 
other by external forces; but if they are HSs, then their 
closest distance of approach is a c below which they face 
an infinite hard-wall barrier. 

The current paper comes as a continuation to the in- 
vestigation of the properties of HC bosons in a CHOCL 
by Sakhel et al. [39(, in continuous space and for a few- 



body system. This is in an attempt to provide further 
manifestation of the condensate properties inside an OL. 
Therefore we ask ourselves various questions: Will an in- 
crease of the HC diameter deplete the SF or BEC inside 
the CHOCL in the case of a few bosons? Can we real- 
ize a MI by increasing the HC diameter to large values 
instead of the OL depth? What is the effect of interfer- 
ence between the condensates in all lattice wells on the 
condensate in one well? We also would like to empha- 
size that the current investigation provides further evi- 
dence for the presence of superfluidity in a 3D CHOCL. 
Previously, Sun et al. [|[ presented striking experimen- 
tal evidence for the presence of SF states in a shallow 
OL. Sun et al. || found that the CF in a deep lattice 
is significantly lower than 1. Further, in a shallow limit, 
the SFF is 100% [3|. Correspondingly, our findings here 
indicate an overall CF of the order of ~ 100% in the 
weakly-interacting regime for N — 8 particles and an OL 
of depth Vq — 10 (in trap units). On the other hand, 
the SFF is ~ 90% in the weakly-interacting regime for 
Vo = 10 and N = 40 particles (see Fig.©. 

The key findings of this paper are as follows: i) the 
most important result is an opposing behavior of the 
global CF and SFF as functions of the lattice spacing. 
Whereas a reduction in the lattice spacing decreases the 
global SFF, it leads at the same time to an increase in 
the CF of the whole system. This is counter-intuitive and 
arises from a distinction between local condensate den- 
sity and global SF density; ii) the SFF is reduced with a 
rise of the HC diameter in correspondence to a decline in 
the single-particle tunneling, earlier reported by Sakhel 
et al.; iii) the interference between the condensates in all 
lattice wells enhances the CF in each well beyond the re- 
sult without interference effects; iv) the energy rises with 
increasing HC diameter whereas the SFF declines. As 
a result, one can conclude that the rise in the energy is 
mostly due to a buildup of the onsite repulsive energy in 
each lattice well which overtakes the drop in the boson 
mobility; v) for a small number of HC bosons N inside 
a CHOCL with a limited number of lattice sites Nl, it 
is possible to achieve a mixed SF-MI state, vi) the local 
and global CF, and SFF, are distinct quantities, vii) the 
principle factor in depleting a BEC in a CHOCL with 
few bosons is the OL. 

The paper is organized as follows. In Sec|Hl we outline 
the methods used for the evaluation of the CF and SFF. 
In Sec lIIIl we present our results and discuss them, and 
in Sec lIVI we present our conclusions. 



II. METHOD 

In this section, we describe the methods used in the 
evaluation of the local and global CF and SFF in each 
lattice well, as well as the global CF and SFF. We con- 
sider N HC bosons confined by a CHOCL of 3 x 3 x 3 
sites. The local CF is computed by separating the one- 
body density matrix (OBDM) of the whole system into 
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its components centered at the various lattice sites. By 
treating the distribution of particles in each lattice well 
as a spherically symmetric cloud (we show justification 
for this) , we diagonalize the OBDM components and hnd 
their corresponding eigenvalues for a particular number 
of sites representative of the whole OL. This is performed 
for two cases, one involving no interference between the 
condensate clouds in all lattice wells, and the other with 
this interference included. The global CF is then ob- 
tained by summing the contributions from individual lat- 
tice wells. The local SFF is evaluated by dividing the 
large cubic OL volume into 27 small cubes of edge d, and 
computing the SFF within the boundaries of each small 
cube. The global SFF is computed as well, and for both 
local and global SFF a winding-number like formula is 
applied. 

For the purpose of evaluating the OBDM and the as- 
sociated eigenvalues, we modified a previously written 
code by Dubois and Glyde [40]. For the evaluation of 
the SFF, we applied the diffusion formula of Pollock and 
Ceperley [4l| . The variational path integral Monte Carlo 
(VPI) method 0, |H| was used to evaluate the SFF of 
the systems treated in Ref.[3!|. We used quantum vari- 
ational Monte Carlo (VMC) to compute the spatial con- 
figurations of the particles in the current systems from 
which the OBDMs are obtained. 



A. Hamiltonian 

The Hamiltonian of the system is given by 

,2 » N 

8=1 1=1 

J^ VintiTi - Tj), (1) 

where Vhoi^i) = \ muj 'ho r 'i lii the harmonic oscillator 
(HO) trapping potential, with r 4 = (x,, J/,, Zj) the posi- 
tion of a particle from the center of this trap, m the mass 
of the particle, and Who the trapping frequency. Upon 
this trap, there is superimposed an OL potential 

Vopti^i) = V [sm 2 (k x Xi) + sin 2 (fc a y 4 ) + sm 2 (k z Zi)] , 

(2) 

with Vo the height of the OL barrier, and fe$ = ir/d 
(i = x, y, or z) is the lattice wave vector with d the 
lattice spacing. For further details refer to Ref.[3!|. The 
interparticle interactions are given by the hard sphere 
(HS) potential 

*wr«) = {~ : r ;;.H , (3) 

where a c is the HC diameter of the bosons, and 
v . — tj \ is the interparticle distance between two bosons 



i and j. a c equals the s-wave scattering length in the low- 
energy and long-wavelength approximation of the two- 
particle scattering problem. 



B. Density matrix and CF at a lattice site 

We begin with the VMC trial wave function of a pre- 
vious publication [3^| given by 

*(W,{R}) = 

N 

J] exp(-arf) ^(r if {R}) J] /(|r, - r,-|), (4) 
where f(\ri — rj\) is the HS Jastrow function [40| 

m = i - (5) 

with a c the HC diameter of the bosons. An important 
condition is that f(r) = if r < a c . Here {r} = 
(ri, r 2 , • • • , r^) is the set of N particle positions, and 
{R} = (Ri, R 2 , • • • , Rat l ) the set of Nl lattice site po- 
sitions, where R„ = (ii + jj + kk)n/d and is abbreviated 
R„ = (ijk). The index n runs from 1 to Nl, where 
Nl is the total number of lattice sites (here 27). Eq.Q 
is then optimized with respect to its parameters [a of 
Eq.dU above, and 8, 7, and a of Eq.® below] as out- 
lined previously [39] . By using the Wannier-like function 
defined by 

N L 

^(r,,{R» = ^^(ri.Rn), (6) 
with </>(r,,R„) given by 

0(r,,R„) = exp[-/3(ri - R„) 2 ] x 

[1 + j(x t - X n f - a{ Xi - X n ) 4 ] x 

\l + l( yi -Y n ) 2 -a{ Vl -Y n f] x 

[1 + 7 (z, - Z n f - a( Zi - Z n ) 4 ] , (7) 

we consider expanding the OBDM into its components 
at each lattice site of position R„. 

To set the stage, we start out by evaluating the total 
density matrix 

p(ri,ri) = J dr 2 dr 3 ---dr N 1 ii*(r 1 ,r2,---,r N ,{R}) x 
^^•••^{R}), (8) 

using the standard Monte Carlo integration approach Q 
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p(ri.ri) = pJ]$ , (r ll r tJj - 1 r EVl {R}) 



c=l 



**(ri,r c2) ---,r cJV ,{R})- 



(9) 



where w c is a configurational weight to be determined 
later, and P is the number of Monte Carlo configurations. 
The subscript c in r C i labels the configuration to which 
particle i belongs. Since according to Eq.© ^'({r}, {R}) 
is a sum over all lattice sites, we can expand /9(ri,r^) 
[Eq.([5|] into a sum of components at positions R„. The 
latter sum p(ri,r[) involves interference components be- 
tween one at R g and all other lattice sites R n7 t g . In this 
study we are concerned with both interfering and nonin- 
terfering ones. On substituting (0} and © into (jSJ), one 
gets 



p(r' 1; ri) = e - ar 'i e - ar i x 



^(^(r^R,, 

,n=l 



'N L 
,n=l 



c=l c Kj l<j 



JV 

n e ~ Qrc2 



0(r CJ ,R„) 

.n=l 



/V 



n /(Ki-rcii) 



(10) 



Let us now choose the configurational weight to be of the 
form 



|*[(r 1 ,r c2 ,.-- J r cW ),{R}]| 2 



= e 



-2arf 



Nr 



N 



^0(ri,R„ 

n=l 

^^(r cl ,R„ 



.n=l 



Kj 



(11) 



which has been chosen in a previous publication [39j. 
Substituting this weight into (fTU|). the exponentials and 
Jastrow terms /(|r G j — r CJ |) cancel out, and one remains 
with 



P(ri,ri) 



— x 



e "'i 

^N L ±<„ -a <\ ^N L 



E^0(ri,R„) E;=i^i,Rn) 



j_ \ - Kj 

Kj 



(12) 



i. Interfering and noninterfering components of the OBDM 

In what follows, we shall separate the density ma- 
trix r' x ) into interfering and noninterfering (isolated- 
islands) components. First, the numerator in the second 
line of Eq. fTl?)) is rewritten 



Nl 



^0(ri,R n ) 



N L 



5>(ri>R«) 
J Lg=i 



E 



(riAJMA) + ^</)(ri,R„)^(ri,Rg) 



(13) 



where the first term in the second line of (| X3[) is the non- 
interfering (n = q) part describing isolated clouds, and 
the second term the interfering part (q ^ n), respectively. 
Next, by using (TT5]l in (|T2"j) . the total density matrix can 
be written as the sum of two components 

N L 

p(r 1: r\) = J2P iQ) (n ~ Rn,r'i - R») + 



N L 



^p (1) (ri - R„,r'i - R n ), 



(14) 



n=l 



where we define the noninterfering density matrix for 
each "isolated" cloud at a lattice site R„ as 



/9 (0) (r 1 -R„,r' 1 -R„) = 



ri,R»)(/)(r / 1 ,R») , 

^-4 r l; r l)> (15) 



Nl 



Lg=l 



and the density matrix involving interference components 
only (n 7^ g) as 



p«( ri - Rn.ri - R„) = 



r , 2 Qr2 ^( ri ,R„) E,^,, X r iiR?) „ , 

1 e 1 — = ./ir' 

^0(n,R 
--3=1 



J(ri,r x ). 



(16) 



For brevity, we have defined the term involving Jastrow 
functions as 



■5 



/(ki 



^ ni^/dri-rc 



(17) 



Note that the order in which the sums of <^(ri,R„) in 
Eg. (|12p are multiplied and divided is very important to 
obtain Eqs. (fT5j) and (fT6|) . First, we evaluated the numer- 
ator (fT3|) by extracting out the interfering (n ^ q) and 
noninterfering terms (n — q). Second, we divided by the 



denominator 



Yet if we considered 



first canceling out the summation Y^ n =i ^( r ij ^-») m the 



numerator of the second line of Eq.(fT2j). with one power 
of the sum in the denominator, we would get a different 
result given by Eg . (|30|) next. 



2. Eigenvalues of the CHOCL-OBDM components 

Assuming now the cloud at each lattice site to be 
spherical, we utilize in what follows the recipe of Dubois 
and Glyde [4(| , which is for the evaluation of the CF in 
an isotropic harmonic trap, to compute the CF in each 
lattice well. Let us first consider the noninterfering case. 
First, one begins by expanding the density matrix (|15[) 
at lattice site R„ into angular momentum components 
as follows: 



Pi (cos j)Pe> (cos 7) sin jd'y 



25, 



21 



1 



(21) 



with 5((> the Kronecker delta function. 

Second, one defines the local natural orbitals (LNOs) 
and diagonalizes the local OBDM to evaluate the eigen- 
values of the OBDM in each well. Hence, similarly to 
Dubois and Glyde [4(| one first redefines the local non- 
interfering density matrix (|18[) using their field-operator 
approach: 



p(°)(n - R n ,v[ - R„) = (#J(n - R„)* (ri - R„)>, 

(22) 

^o( r i — Rn) being a local field operator which annihilates 
a particle at position ri — R„. Then one expands $o( r i — 
R„) into a set of local single particle states Xi( r i — Rn) 
around R„: 



* (ri-R„) = ^x*(ri-Rn) 



(23) 



a\ n ^ being a bosonic annihilation operator acting inside 
well n at position R„. The \i are taken to be orthonor- 
mal such that 



^'(n-R^-RO = 



E 

1=0 



( 2i + 1 ) o (M \ (0) 



-in 



ft(^)pf (|r 1 -R n |,|r' 1 -R n |), 



(18) 



J dr lX *{ri - R„)xj(ri - R 



(24) 



Subsituting (f2"3"]l into (j2"2")l above, and using the usual 



condition (a^^a^) — N^Sij, one gets 



where P e (x) is the Legendre polynomial of order £ and p (o) (ri _ Rn; r ;_ Rn) = £ x*(ri-R«)Xi(ri-R„)ivf 



(„) _ (n - R n ) (r[ - R M ) 

- In-Rnl ' |ri-Rn| ! 



(19) 



the cosine of the angle between (n — R„) and (r^ — R„). 
The pf^ corresponds to the angular momentum compo- 
nent of p(°^ and is obtained by multiplying both sides of 
(jl8]l by Pi (it,! r ' ) an d integrating over the solid-angle ele- 

(n) 

ment dili = 2ir sin 7^7, considering that cos 7 = u ri T , . 
This yields 



pf> (|r 1 -R n |,|r' 1 -R n |) = 

J dKliPt p(°)(n - R„, ri - R„), 



(20) 



where the orthogonality condition for Legendre polyno- 
mials [45[ has been applied: 



(25) 

Multiplying both sides of (PZBT) by Xj( r i — R-n) from the 
left and x*( r 'i — R-n) from the right, then integrating over 
Yi and r[, this yields 

J dndrixj-(ri - R„)p (0) (ri - R^.r'i - R„) x 
Xjfri-Rn) = NiS tj . (26) 
The eigenvectors x«( r i ~ are the LNOs with eigen- 

(n) 

values N> '. Making now use of the spherical symmetry 
of the system, one further defines 



Xi(rj - R„) = ^(|ri - R n \)Y tm (e, <p), (27) 

i = (q£m) being a state, £^(l r i — Rn|) being the local ra- 
dial wave function, and Yg m (6, <p) the spherical harmonic 
function with 6 and ip defining the angular position of 
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rj — R„. It must be noted, that in we drop the 

quantum number m so that N- n ' = N^' because of az- 
imuthal symmetry in each well. Substituting (|27[) and 
(fT8"|) into (f2"5l) . using the addition theorem for Legendre 
polynomials [45j 



P e (cosj) 



4tt 



m=+£ 



21- 



T £ r, m (0^)y/ ro (0V) (28) 



cos 7 = sin sin 0' cos(<^ — y?') + cos0cos#', 



with 



and the orthogonality condition (f21j) . this eventually 
yields the expansion of the angular momentum compo- 
nent pp|) into the radial functions £ q e 



,(°) 
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(Iri-Rnl, 



ri-Rnl) 



Rn|)JV, 



(n) 



(29) 



where AQ™^ are the eigenvalues of the natural orbitals 
Xi( r i — R-n) at lattice site n and positions R„ = (ijk). 
As such, the natural orbital Xo( r i — R-n) whose eigenvalue 
is Nqq\ is defined as the condensate orbital. Nqq/N is 
then the CF at site n and position R„ with respect to 
the overall number of particles N. One then proceeds 
similarly to Dubois and Glyde in the diagonalization 

(n) 

of the OBDM in order to compute the eigenvalues N^ e , 
and consequently the CF, except that we do this here for 
one lattice site at a time. 

By considering the additional interference between a 
condensate at lattice site R n and all-neighbor lattice sites 
R q 7^ R„, the CF at lattice site R„ is computed using the 
same previous procedure, except that we use all terms of 
the density matrix (|12[) which contains the interference 
between the condensates at the various lattice sites, in- 
stead of the noninterfering p(°> only. That is, one writes 



-Rn.ri 
(ri,R„) 



Rr, 



(30) 

where p = p (0) + Consequently, Eqs.([THI), (f20]>. and 
(|29|) are simply rewritten by replacing p^> with p in all 
terms. Eventually, one gets instead of (|29|) 



C. Definitions of CF in our systems 



To set the stage, we define nujjA — Nqq /N to be the 
CF in the lattice well centered at R„ = (ijk) with respect 
to the total number of particles N, and without all- well 

interference, and 



Hijk) 



= Nqq/N as the corresponding 
CF with all-neighbor interference, labelled by a star in 
each upcoming figure legend. The nujiA and fi(ijk) are 
computed as explained in Sec lIIBI 

Hence n( oo) is the fraction of particles in the lowest 
orbital of the central cell, n(oio) that of cell (010), etc. 
In fact, the natural orbital is spread out over the entire 
lattice and is not localized in one lattice site. That is, one 
has a lowest orbital at all lattice sites (ijk) if they display 
a CF. This spreading is a Wannier sum over individual- 
cell natural orbitals. Since the CF at site (ijk) is with 
respect to the total number of particles N, the total CF of 
the system in the 3x3x3 cubic lattice can be computed 
using 

ijk 

"•(ooo) + 6rt( io) + 12n(on) + 8n (111 ). (32) 

Note that n(oio) is ideally identical to ri(ooi)) "(100)1 
«(-ioo), etc. Similarly 7i( n) is equal to n (10 i), "(no), 
n(_no), etc; the same holds for Tinny The large cube 
divided into 3x3x3 smaller equal-sized cubes has 6 
lattice sites at its face-centers, 8 at the corners, and 12 
at the bisecting points of its edges. 

The CF in each lattice well is then explored as a func- 
tion of a c for the three lattice wave vectors k — tt, 1.2tt, 
and 1.47T. The latter results are compared with the cor- 
responding CF flUj]A . 



D. Superfluid fraction 

The global SFF is computed using the variational 
path integral Monte Carlo (VPI) technique I42|. [431 via 
the diffusion formula of Pollock and Ceperely [4l[ : 



El 
p 



Dp 
A)' 



(33) 



where Do 
stant and 



Ti 2 /(2m) is the "quantum diffusion" con- 



Pt (\T! -RnUri -Rnl) = 

^2UQri - R»|)^(|ri - Rn\)N$\ (31) 



where N^' and ^ are similar to N$ and but 
including contributions from the interference with all- 
neighbor lattice sites. 



2d/3N 



N 



(34) 



with d the dimensionality of the system, N the number 
of particles, imaginary time j3 ■<-> l/(ksT) with fcs the 
Boltzmann constant, and T the temperature. Here (• • •) 
denotes a configurational Monte Carlo average, is the 
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initial position of a particle i, and YMi the destination 
of the particle after a "time" (3, where M is the num- 
ber of VPI time slices. In this regard, we evaluate the 
configuational average D p over a number P of VPI con- 
figurations: 



D 



1 1 



p,VPI 



2dfiN P ^ 



N 



^(r ci - v M ci) 



(35) 



For d = 3, and using units of the trap a^ = yft/ (muiho) 
and ?iw; lo for length and energy, respectively, the SFF is 
then recast into the form 



1 1 p 
Ps_ _ 1 1 

p ~ 3MtNP^ 

c—l 



N 



^(r ci - r Mcl ) 



i=i 



(36) 



where r — ► r/a^o, P — > /3Tiu>ho, and t — > pTvUho/M — 
TiLOhoj '(MksT) is the "time step" with M the number 
of time slices. We would also like to remind the reader 
that we are aiming at presenting qualitative rather than 
quantitative results. The parameters used were M — 
120, r = 5 x 1(T 3 , and N = 40. The number of MC 
blocks was ~ O(10 3 ), similarly for the number of MC 
steps for each block. 

In order to compute the local SFFs in the individual 
CHOCL wells using Eq.(|3"6"]l. we divided the large vol- 
ume of our cubic OL into 27 smaller cubic cells of equal 
sizes d 3 . The cube edge was equal to the lattice spacing 
d = ir/k. The central cell was centered at the origin 
(000) with its faces parallel to the coordinate planes, the 
rest of the cells being centered at the other lattice sites. 
By restricting Eq. ([36]) to the boundaries of each cell vol- 
ume, we were able to compute the SFF (p s /p)ujk) in each 
lattice well. For each cell (ijk), we compute (ps/p) a (ijk) 
using only the particles which are positioned inside the 
cell according to the boundary conditions 



1\ 7T / 1\ 7T 

-2)k <X< \ l+ 2)V 

1\ 7T / 1\ 7T 

. 1 \ n (i 1 \ n 

k -~2 k <X< ( k+ 2 T 



That is, considering (f3"T|) the local SFF is 



(37) 



E. Optical densities 

In this article, we are also concerned with the average 
VPI integrated optical 2D density, (n 2 D(x, y)). The in- 
tegration is of the total density along the z— axis. That 
is, we first define a density ri2D,c(x,y) for each VPI con- 
figuration c such that 



n2D,c(x,y) = 



|*c(z, V, z)\ 2 dz, 



(39) 



where ^ c (x,y, z) is the wavefunction [Eq.@] obtained 
for a certain VPI configuration c. Then one takes the 
average over all P configurations 



(ri2D(x,y)) 



-t 

P ^ 



n 2 D,c( x 



(40) 



Next to this, we also display at some point the ID optical 
density, (nifl(i)), which is obtained by integrating the 
2D density along one of the axes: 



(n 1D (x)) 



+ 00 



dy(n 2 D(x,y)). 



(41) 



We therefore would like to point out, that in Ref.[39( 
we mistakenly wrote that tiid{x) — (fi2D(x,y = 0)), 
whereas it should be as in Eq. (|4T|) . Nevertheless, in the 
main text of Ref.(39| and in its figure captions, it was 
clearly indicated that ti\d{x) is the integrated ID optical 
density, i.e., as in Eq. (|4T]) above. 



F. Units 

As in Ref.[39| we use units of the trap, ai lo = 
\Jfij (muJho) and Tuouo for lengths and energies, re- 
spectively, all throughout our calculations. That is, 
r ->■ r/a ho , V -> V /(hu) ho ), a c -)• a c /a ho , V ho {r) ->■ 
Vho(r)/(hu h o) = (l/2)r 2 , k -)• ka ho . The optical den- 
sities {ri2D{x,y)) and {nijj{x)) are in units of and 



respectively. 



III. RESULTS AND DISCUSSION 



\PJ{ l3 k) \PJ \ N ^jk)) 

In this paper, {p s / p)(ijk) is computed using the latter 
equation, but sometimes it is found more reasonable to 
compute it with respect to the total N, particularly if 
(N)njk) — > in the local MI phase. That is, we use 
(f31)f plus (f3"T|) . but without multiplying by the factor 
N/{N (m ). 



In this section, we present the results of our numeri- 
cal calculations. The same systems are considered as in 
Ref. [39(| . The latter are V HC bosons confined inside a 
cubic OL of Vl =3x3x3 lattice sites plus an exter- 
nal harmonic trap. For further details about the system 
and the computational setup, please refer to the previous 
publication |39j . 

We explore BEC and superfluidity both from a global 
and a local perspective. The local perspective, pertains 
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FIG. 1: Density map for a system of N = 8 bosons, Vo = 10, 
and k = l.On in the same confining geometry as Fig. [3] below. 
This map is obtained by integrating the three-dimensional 
density along the z— axis. 



to their fractions in each lattice well (or cell), and the 
global one pertains to these fractions for the whole system 
(all lattice wells). 

From the local perspective, we begin with a brief in- 
vestigation of the ideal Bose gas, in which we show that 
the lattice wells have a low CF for N — 8 particles when 
viewed with respect to the global N. Next, we demon- 
strate the dependence of the CFs nujf.) and fi(yfe) on the 
HC diameter a c , where an important role for the interfer- 
ence between the lattice wells is revealed. Similarly, the 
global and local SFFs, p s /p and (p s / p)ujk), respectively, 
are explored as functions of a c . 

From a global perspective, we chiefly find a counter- 
intuitive, opposing behavior for the global CF and SFF 
as functions of k — w/d. Peculiarly, the SFF drops with 
increasing k while the CF rises. Additionally, we seek a 
MI state by going to larger interactions or OL depth, and 
show that our systems display a coexistence of SF and 
vacuum MI regimes. 

Finally, we also present results for the momentum dis- 
tributions for some of our systems. The latter reveal 
a diffractive structure which persists deep into the MI 
regime [46| suggesting that superfludity [46[ in a CHOCL 
is found even in the extremely repulsive regime. 



A. Numerics 

For the CFs at the lattice sites (ijk), the configura- 
tions of previous VMC simulations are used [39[ . The lo- 
cal CF is computed as in Sec lIIBl whereas the global one 
by the additional Eg .([32]). The global SFF is computed 
using Eq. (j5rJl) and the local SFF by an additional appli- 
cation of the boundary conditions (|57|) . We use N = 8 
particles with V = 10, and N = 40 with V = 10 and 
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logio(a c ) 



FIG. 2: CF n(jjfc) = N^^/N at the central cell (open cir- 
cles) and (010) (solid circles) as a function of a c in the ex- 
tremely dilute regime. The system is a HC Bose gas of N = 8 
particles, Vo = 10, and k = 1.07T in the same confining geome- 
try as Fig. [3] below. A log-log scale is used in order to provide 
a clearer view of the details. The o c , k and Vo are in units of 
dho, a^l ancl hwho, respectively, where a ho = ^/h/(mwho)- 



once with Vo = 20. The total number of VPI time slices 
used is M = 120 with a time step r - O(10" 3 ). The 
number of MC steps and MC blocks is - O(10 3 ). 

In order to justify the method outlined in SecJlTl ne- 
cessiating the presence of spherically symmetric clouds 
centered at each lattice site, we present in Fig.Q] a 
map for the integrated VMC density along the z— axis 
(ri2D(x,y)) for the system of Fig. [3] below. One can see 
clearly that the distributions of the atoms about the lat- 
tice sites are spherical at the edges of the trap. Thus it 
can be safely stated, that for a few bosons spherical sym- 
metry is displayed by the lattice clouds, in spite of the 
presence of an external harmonic trap. Yet by increasing 
the number of particles, it is anticipated that the clouds 
at the edges of the trap will be deformed from their spher- 
ically symmetric form, as the atoms are pushed away 
from the trap center towards the edges of the trap by 
their mutual HC repulsion. This is especially the case, 
when a c becomes large. Consequently, in the absence of 
spherical symmetry, the method of Sec|TT] cannot be used 
anymore, except for the clouds which may nevertheless 
preserve their spherical form. If one uses a large CHOCL 
with, say, Nl = 5x5x5 lattice sites, one must consider 
a weaker harmonic trap than the one used here in order 
to maintain spherically symmetric clouds at the edges of 
the trap. 



B. Condensate fraction and interference effects 

In this section, we explore the local and global BEC 
in our CHOCL systems. Only five lattice sites, (000), 
(001), (010), (011), and (111), which present the whole 



9 



CO 

o 



o 

. - — I 

o 

a3 
cfcl 
O 

d 

o 

o 
O 



0.16 
0.14 
0.12 
0.1 
0.08 
0.06 h 
0.04 



0.14 
0.12 
0.1 

0.08 - 
0.06 - 
0.04 



0.14 
0.12 
0.1 
0.08 
0.06 
0.04 



k = I.Ott h— h 

I.OtT* t-B~i 

R„ = (000) 
AT = 8 
V o = 10 

l^H— B--E „ 



k = 1.27T 

1.2tT* :~V- 



le- 



ft V- 

-e-a v -v 



k = 1.4tt^ 
1.4-7T* ♦ 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 

a c (units of a^ ) 



FIG. 3: CFs ri(ooo) and n(ooo) at lattice site R n = (000) versus 
the HC diameter a c for systems with and without interference 
effects, respectively, and at the various k— values indicated. 
The systems are the same as in Ref . [39l] : N = 8 HC bosons 
confined in a cubic OL of Nl = 3x3x3 lattice sites plus 
external harmonic confinement. The height of the OL barrier 
is Vo = 10. The fifooo) values are labelled by (*) in each 
legend. Solid squares: k = I.Ott; open squares k = I.Ott*; 
open circles: 1.2tt; open triangles 1.2tt*; solid triangles: 1.4tt; 
diamonds 1.4tt*. The k and Vo are in units of o^ 1 and Tuoho, 
where an = \f%} (moj ho ). 



CHOCL system are considered. The role of interference 
between the lattice wells is revealed by a comparison be- 
tween nujfA and h/ijk). 



Local CF in i 



BEC 



We begin with an investigation of the ideal Bose gas 
for a c values ranging from 10 -16 to 10 -3 and show that in 
this regime the local CF is almost stable. The a c values 
range from an extremely dilute to dilute regime. Fig. [2] 
displays the local CF at the central cell (open red circles) 
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FIG. 4: Same systems as in Fig- El but without interference 
data. The k and Vo are in units of a^ o and Kuiho, respectively. 



and (010) (solid blue circles). We get n( 000 ) ~ 15% and 
"(oio) ~ 5.5% of the total N, which in themselves are 
small CFs. It should be emphasized that the optimized 
VMC parameters of the trial wave function Q are ex- 
actly the same for all a c in the dilute regime. This is 
because in the dilute regime, the wave function of the 
system [Eq.Q] is essentially a product of the Wannier 
sums Eq.© with little effect coming from the Jastrow 
function. 



2. Local CF in the interacting BEC 

Fig. [3] displays the CFs n( oo) and ri(ooo) (starry la- 
bels) in the central well R„ = (000) versus a c at the 
indicated values of k. From the top frame to the bottom 
frame we display for n(ooo) : (solid squares) k = I.Ott; 
(open circles) 1.2tt; (solid up triangles) 1.4tt. The sys- 
tems are the same as in Ref. [39j| . For n( oo) we display: 
(open squares) k = I.Ott*, (open down triangles) 1.2tt*, 
and (diamonds) 1.4tt*. Fig. [4] is the same as Fig.[3l but 
without interference data. The goal of Fig. |4] is to give 
further manifestation to the dependence of n(ooo) on k. 
One can see that for all values of k in the latter figures, 
n,(ooo) an d ^(ooo) decrease with an increase in a c . Further, 
the n(ooo) lies higher than ri(ooo) for all k, indicating that 
interference between the lattice wells boosts the CF in 
each well. In Fig.|H the values of n( 000 ) seem to rise with 
increasing k in the range 0.05 < a c < 0.25, although for 
k = 1.2-k and 1.4tt the values are very close to each other. 

Fig.[5]is the same as Fig. H but for R„ = (001), (010), 
(011), and (HI), ordered respectively from the top frame 
to the bottom frame. The same labels are used as in Fig. 
13] However, whereas n and h decrease with increasing a c 
for (001) and (010), they surprisingly increase for (111) 
and (011). Some condensate has tunneled away from the 
sites R„ = (000), (001), and (010) towards (011) and 
(111) as a result of an increase in the mutual repulsion 
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FIG. 5: As in Fig. [3] but from top to botton for lattice sites 
R„ = (001), (010), (Oil), and (111). The k and V are in 
units of a^ and Tiu>ho, where aho = v^7 (mcoho)- 



of the bosons which drives them away from the center of 
the trap. A similar expulsion of the BEC towards the 
edges of the trap has also been encountered in the simple 
harmonic trap without an OL 5?} . This phenomenon 
corresponds also to a rise in the occupancy of the lattice 
wells at the corners of the CHOCL with increasing a c , 
as demonstrated in Sec lIII El later on. Otherwise, Fig.[S] 
displays the same features as Fig. [3] regarding the values 



of k and the interference effects. The rise of 
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increasing k is even more pronounced in Fig. [5] One can 
explain this by approximating each lattice well by a HO 
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FIG. 6: CF ri(ooo) versus the HC diameter a c for a system 
of N — 40, Vo = 20, and k = 1.0-7T in the same trapping 
geometry of Fig. [3] The k and Vo are in units of o^ 1 and 
TiuJho, respectively, where a^o = v/S/(tHW(,„). 



trap, whose trapping frequency ujk rises with increasing 
lattice wave vector k = w/d. As a result of increasing uik, 
the level spacing ~ TiuJk becomes larger, thus increasing 
the energy needed to excite atoms out of the ground state 
HO level. Consequently, more particles find themselves 
forced to condense into the lowest "harmonic oscillator" 
level. (We did not include in Fig.[5]n for k — I.Ott* in 
order not to clutter the figures.) 

Another explaination is as follows. Going back to Fig. 
15 of Ref . (39| . one can see that the occupancy of the cen- 
tral well (000) and the first neighbor (00 —1) decreases 
with increasing k. [Note that (00 —1) belongs to the same 
lattice site family of (010) and (001).] Consequently, the 
reduction in the occupancy reduces the density in these 
lattice wells, and henceforth the onsite interactions be- 
tween the bosons. This in turn reduces the depletion 
effect of the HCs causing a rise in the CFs with increas- 
ing k as observed in Figs. [3j— [5j It is further noted that 
n(ooi) is changing with a c at a lower rate than n( oo) • E.g, 
n.(ooi) f° r k — 1.07T (solid squares) in Fig. [S] declines from 
~ 0.055 down to ~ 0.04 over the whole range of a c in the 
figure, whereas n( oo) for the same k (see Fig. [4]) decreases 
from ~ 0.12 to ~ 0.04. The same feature is observed for 
R„ = (010) when compared to (000). In contrast, nmi) 
and Ti(oii) increase by ~ 0.01 in Fig.[SJ Hence, the rate 
of change of the ^(yfe) with a c , whether Ji(»jfe) is rising or 
declining, becomes smaller as the lattice well approaches 
the edges of the trap. This is attributed to the fact that 
the energy cost required for the bosons to tunnel through 
(or between) two neighboring CHOCL wells increases as 
the lattice wells approach the edges of the background 
harmonic trap. This is because each lattice well is super- 
imposed on the harmonic trap. Therefore the potential 
energy minimum of each lattice well at position R n from 
the trap center, is equal to the HO energy (l/2)|R n | 2 
and rises with increasing |R„|. In addition, we present 
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in Fig.[S]the behavior of the CF ri(ooo) f° r N = 40 parti- 
cles, Vo — 20, and k = n. Again, the BEC gets depleted 
with increasing a c , except that the CFs are substantially 
higher than for N = 8. 

It must be emphasized, that the values of the CFs in 
Figs. |3H5] fluctuate with a c due to the small number of 
particles, N = 8, used which causes real physical fluctua- 
tions to be large [48j]. Nevertheless, the trend in the CFs 
clearly indicates depletion with increasing a c for some 
lattice cells and a condensate buildup in other cells. It is 
further noted, that in general the CFs for N — 8 are sig- 
nificantly lower than 1, bringing this result in line with 
the earlier findings of Sun et al. Q. For N — 40, the 
data varies smoothly with a c and the statistics are good, 
demonstrating that a higher N leads to lower condensate 
number fluctuations inside the OL. 



3. Global CF n 



regime and remain in this state even up to very large HC 
repulsion. Further, the global CF is boosted with the rise 
of k (the local-confinement strength ujk oc k = n/d) in 
each CHOCL well. The OL introduces a local depletion 
effect at lower k < lAn when there is almost no role for 
the HC repulsion in the global depletion. An important 
finding, then, is that the BEC in a CHOCL with a few 
bosons is chiefly depleted by the OL. Therefore, whereas 
the global CF no remains constant as a c is changed, the 
local CF is redistributed at the various lattice sites. In 
contrast, Ramanan et al. [33[ found that for a ID Bose 
gas trapped by an OL plus a weak external harmonic 
trap, the total CF is depleted substantially with an in- 
crease of the onsite repulsive energy U. Using the BHM, 
van Oosten et al. [12] demonstrated that the CF in a 
2D and 3D OL decreases with increasing parameter U /t, 
where U is the contact interactions strength and t the 
hopping amplitude. 



Using the previous results of Figs. [3]- [5] for the non- 
interfering case only, we can estimate the global CF no 
using Eq. (|3"2"|) . The goal is to distinguish between conden- 
sate local and global behavior. Fig. [7] displays no for the 
systems in the latter figures at k = I.Ott (solid squares), 
1.27T (open circles), and 1.47T (solid triangles). It is noted, 
that no is practically constant and not sensitive to the 
changes in the repulsive forces via a c . This is a peculiar 
result, suggesting that the CF is globally conserved. In 
contrast, the decline of ^(yfc) with increasing a c at the 
lattice sites (000), (001), and (010) in Figs. H- M is a lo- 
cal effect and was compensated by a corresponding rise at 
(011) and (111). The n for k = I.Ott is ~ 80% for all val- 
ues of a c , whereas for k = 1.2n and 1.47T, no ~ 85% and 
~ 100%, respectively, for all a c . The values of no are high 
because for N — 8 particles the systems are in the dilute 
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FIG. 7: Global CF no as estimated by Eq.([22J for the systems 
in Fig.[3]at k — I.Ott (solid squares); k = 1.2tt (open circles), 
and k = 1.47T (solid triangles). The k and Vo are in units of 
a^l and huiho, respectively, where an = \/h/(mojho) 



C. Superfluid fraction 

In this section, we explore the SFF in our CHOCL 
systems both as a global and as a local quantity in in- 
dividual lattice wells. We consider systems with N = 8 
and N = 40 bosons in the same confining geometry as in 
Fig. [3] with the same k— values and Vq — 10. The SFF is 
computed globally using Eq. (|3"o) . and locally for a c < 0.3 
by Eq. (|38[) and an additional application of the boundary 
conditions (f37|) . For the larger a c ~ 0(1), the SFF was 
computed with respect to the total N as outlined later 
below. We did not compute the SFF for the systems 
with N = 8 particles for the range of a c < 0.3 consid- 
ered so far because of the high statistical fluctuations in 
the values of the VPI MC block averages. However, at 
extreme values of a c beyond 1.00, these fluctuations are 
much lower, and therefore it was reasonable to report the 
SFF for N = 8 as in SecUHE] below. 



1. Global SFF p s /p 

Fig. [5] displays the global SFF, p s /p, versus a c for 
the latter systems at the indicated values of k. The la- 
belling is as in Fig.[7J One can see that the global SF 
is depleted with a rise in a c . Further, the values of p s /p 
decrease with increasing k, and one can conclude that a 
higher k (lower lattice spacing) enhances the depletion 
of the global SF in an OL. This is our chief result and 
is counter-intuitive to the rise of the global CF with in- 
creasing k as in Fig. [7] Retrospectively, this could also 
be explained as before [3^ |: a smaller lattice spacing in- 
creases the localization of the bosons inside each lattice 
cell and causes a reduction in their superflow. A reduc- 
tion in superflow means a smaller number of particles 
possessing enough kinetic energy for tunneling from one 
well to another. Further note that remarkably even at 
very strong repulsion between the bosons, superfluidity 
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FIG. 8: Global SFF p s /p versus the HC diameter a c as com- 
puted by the diffusion formula of Pollock and Ceperely [4l|] 
given here by Eq. (|36[l . The system considered is TV = 40 HC 
bosons in the same trapping geometry of Fig(3] and the indi- 
cated values of k. The same labels are used as in Fig. [5] The 
k and Vo are in units of aZ a and hu^o, respectively, where 
aho = ^jTi/{mu ho ). 



is still present (p s /p ~ 40% for a c = 0.3.) In contrast 
to the global BEC in Fig. [3 a c plays a chief role in the 
global depletion of the SF in Fig. |SJ 



2. Local SFF (p s /p)(ijh) 

As for the SFF in each CHOCL cell, {p s /p)(ijk)i it re- 
veals a different behavior as compared to its global char- 
acter in Fig. IU Fig. [9] displays (p s /p)(ijk) vs - °c for the 
same systems of Fig. [FJ where for the top frame to the 
bottom frame: k = I.Ott, 1.2n, and 1.47T, respectively. 
Solid squares: (ijk) = (000); open circles: (010); solid 
circles: (011); open triangles: (111). Inside each cell, 
the local {p s /p)(ijk) is computed using Eg. (I3"5)) plus the 
condition (f3"T)l . 

The same features are observed in all three frames: 
(p s /p) (on) an d (Ps/p)(m) rise with increasing a c . The 
(p s / p) (oio) displays an initial weak decline, but then 
it stabilizes somewhat after a c = 0.2. For the center 
(000) the decline is more pronounced than in (010) up 
to a c = 0.16, after which it stabilizes somewhat. In fact, 
the rise of (p s /p)(on) an d (p s /p)(in) with increasing a c 



conforms to the rise of the CFs 



Hon) 



and 



Hill) 



with 



in Fig. [5] (although for a different number of particles), 
whereas the decline in (000) corresponds to that of ri(ooo) 
in Fig. SJ It can therefore be argued, that locally the be- 
havior of the SF is isomorphic to that of the condensate. 
Although N is different, it is important to emphasize 
that the local SFFs in Fig. [9] are of the same order of 
magnitude as the local CFs in Fig.[5j The reduction in 
{Ps/p) (ooo) i s m une with the previous finding [39|, that 
a rise in a c reduces the single-particle tunneling ampli- 
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FIG. 9: As in Fig. El but for the local SFF (p«/p)(ij*0 at 
the indicated values of k. The cell boundaries are given by 
Eq. (PT)) . The (p s /p)(ijk) is computed via Eg. ((55)) by using the 
number of particles Nuj^) inside each cell instead of the total 
number of particles N. Solid squares: (ijk) = (000); open 
cirlces: (010); solid circles: (011); open triangles: (111). The 
k and Vo are in units of a^ 1 and Tuoho, respectively, where 
aho = \/h/(muj ho ). 



tude J between, e.g., the central (000) and first neigh- 
bor (010) cells. The superfiow is therefore suppressed 
between those lattice wells because of an increased local- 
ization at the larger a c . 

By inspecting Fig. HI one notes again a similar 
counter-intuitive feature as for the global values; that for 
each cell the local SFF drops with increasing k contrary 
to the corresponding local CF in Fig. [5] Another pecu- 
liar feature in Fig. |H1 is that all four curves in each frame 
intersect almost at the same value of a r where there is 
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0.35 



FIG. 10: Average VPI energy per particle, (Eypi/N), versus 
the HC diameter a c for the same systems of Fig. [8] The k 
is in units of a^ 1 where aho = \/Ti/ (mojho). The (Eypi/N) 
and Vo are in units of hojho- 



almost an equal distribution of the SF at all lattice sites. 
The SF migrates from the cells near the center of the trap 
towards the corners of the CHOCL with the increase of 
a c . This feature is also elaborated in Sec lIIIEl below, 
where systems with extremely repulsive HCs are investi- 
gated. 

The effects of an external harmonic trap are as such 
to introduce a new behavior for the superflow of bosons 
between the wells of an OL. As the repulsive interactions 
become stronger, the SF is "expelled" towards the edges 
of the trap analogous to the condensate in a simple har- 
monic trap [4fj, |47|. However, in a CHOCL parts of the 
condensate still remain situated in the lattice wells close 
to the center of the trap. 



D. Energies 

Fig. [10] displays the behavior of the average VPI en- 
ergy per particle, (Eypi)/N, as a function of the HC 
diameter a c for the same systems in Fig- [8] Whereas 
(Eypi/N) rises with increasing a c , p s /p decreases (see 
Fig- EJ ■ The rise in the energy is chiefly due to the rise in 
the average onsite interaction energy {U( pqr }) with (Xq , as 
found earlier in Ref.[3!|. Therefore, the later rise over- 
whelms the drop in kinetic energy, i.e. mobility, in favor 
of an increase in boson localization. 



E. Can we get a pure Mott insulator? 

In this section, we show that for a few bosons in 
a CHOCL with a limited number of lattice sites, MI 
regimes can be achieved by increasing the HC diameter 
of the bosons to large values. It is found that only a part 
of the system will be in the MI state where a number 



of CHOCL wells display the absence of superfluidity. We 
show that our systems are different than the ones usually 
treated by the BHM for a much larger number of par- 
ticles N. For example, Hen and Rigol [l3| explored the 
phase diagram of a HC BHM on a checkerboard superlat- 
tice. According to Hen and Rigol [UHlij], bosons inside 
a CHOCL display coexisting SF and MI phases. Mixed 
SF-MI regimes have also been experimentally reported 
by Spielman et al. [TH, H|| in 2D atomic gases confined 
in harmonic OL potentials. In the work of Jaksch et al. 
[49| . a checkerboard SF-MI phase was observed in an OL. 
Their most important finding was that a BHM could be 
realized by the dynamics of bosons in an OL. But this was 
found for a much larger number of particles and lattice 
sites than ours, using a discrete space approach, whereas 
our investigation is restricted to a few bosons and lattice 
sites in continuous space. 

Hen and Rigol [Til ] explored the ground state prop- 
erties of HC bosons in 2D and 3D OLs confined by an 
external harmonic trap. It was demonstrated that a MI 
state usually displays a flat density along several lattice 
sites, as in their Figs. 7 and 14, where the density is dis- 
played as a function of the distance from the center of 
the trap. However, a limited number of lattice sites as 
in the present work here does not furnish the ground for 
obtaining a flat density profile. 

In the investigations of Hen and Rigol EE!, Si] and 
Jaksch et al. [4!| the interactions are described by contact 
potentials and have no range (a c ) as in our case. Thus, 
the BHM cannot be used here to make predictions about 
a possible SF-MI transition. In fact, a large increase 
in a c causes —as shown next— the bosons to be driven 
out of the center of the trap towards the corners of the 
CHOCL, instead of making them prefer to occupy sites 
individually as in the BHM. This is because when the 
volume of each HS boson is increased, they have no other 
way but to increase their minimum separation as they are 
not able to approach each other by a distance less than 
a c , prescribed by the Jastrow function (JSJ. Therefore, 
a homogeneous distribution cannot be achieved. Hence, 
increasing a c largely in our systems does not yield a pure 
MI state throughout the whole lattice wells. We are thus 
very much inclined to say, given the above information, 
that our systems will retain an SF even in the strongly 
interacting regime. 

In the BHM, on the other hand, a large increase in the 
repulsive contact potential makes it energetically more 
costly to hop from one lattice well to another, but it 
does not increase the volume of the HS bosons. A ho- 
mogeneous distribution is thus naturally possible in the 
BHM as the bosons seek the lowest energetic configura- 
tion in the OL. 

Further, in contrast to Hen and Rigol fl5j ] . the vacuum 
surrounding our CHOCL does not contain empty lattice 
sites, and therefore cannot be associated with an empty 
MI state. The empty MI regime was mentioned earlier 
fl~n| . However, inside the CHOCL it will be demonstrated 
that vacuum MI states are possible. 
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VPI, N = 8,V = 60 




FIG. 11: Integrated VPI density along the z— axis (ri2D(x,y)) 
[Eq.gO])] of a system with TV = 8, a c = 0.02, k = l.Ovr. The 
trapping geometry is the same as that of the systems in Fig. 
except that the OL potential barrier is Vo = 60. The a c , k 
and Vo are in units of aho, cl^I and hu>ho, respectively, where 
a,ho = v^7 { mL °ho)- The x and y are in units of a^ D , and 
(ri2D(x,y)) is in units of a~^. 

In order to check for the presence of a SF phase in 
some parts of our systems for large a c , one could com- 
pute U/(zJ) between two neighboring lattice sites. This 
should help us shed more light on the state of the cur- 
rent CHOCL systems, because it is hard to distinguish 
between a SF and MI by only looking at the densities. 
The identification of MI domains only by means of the 
density was found to be inaccurate, as outlined earlier by 
Rigol et al. [l4[ (and references therein). Earlier, it was 
found [39| that for a c up to 0.3 our systems were still in 
the SF phase and that U /(zJ) began to stabilize for large 
a c . For the tunneling from the central cell (000) to the 
hrst nearest neighbor (00 —1), our calculations yielded 
a U/(zJ) < 5.814 (see Fig.6 in Ref.fli]). However, for 
the much larger a c values used next, an evaluation of 
U I (zJ) was not possible since a VMC ground-state wave 
function was not obtainable anymore for a c > 1. This is 
because VMC reweighting [40| breaks down at these large 
interactions. Instead, we resorted to evaluate the VPI SF 
fraction in each lattice well as was done in Sec lIII CI 

1. Large Vo 

First, we tried to obtain a pure MI by an increase of 
Vo- In Fig. [TTJ we show such a case of high Vq. There, we 
display the integrated VPI density (ri2D(x,y)) [Eq. (l4*Dl ] 
for a system of N = 8, Vq = 60, a c — 0.02, and k = l.On 
in the same trapping geometry as that of the systems 
of Fig. [31 The corresponding Fig. [T^] displays the inte- 
grated ID VPI density [{n 1D (x)} (open circles), Eg. ([41])] 
of Fig. [TT] along the x— axis. The solid line is the inte- 
grated ID VMC density obtained by Eq.flUJ as for VPI. 
It is a VMC best "fit" to the VPI data, obtained by man- 
ually fitting the parameters of the 3D VMC trial func- 
tion [Eq.flD] to the ID VPI density. The goal was to 
use the "fitted" VMC trial function to compute the ra- 
tio (C/(ooo))/[6(J(ooo)^(ioo))] [39] between the central cell 




x (units of dh ) 



FIG. 12: Integrated VPI ID optical density (n 1D (x)) [Eg. ([TO] 
of Fig. [11] (open circles) along the x— axis. The solid red line 
is a VMC best "fit" to the VPI data using the trial wave 
function Eq.Q. The inset is the same figure; but taken for 
a much smaller density range in order to reveal the overlap 
of the wave function between the lattice sites. The a c , k 
and Vo are in units of aho, and huJho, respectively, where 
a ho = ^Jh/(mu) ho ). 



and a first neighbor cell in order to make a decision about 
the state of the system. It was found that the latter ra- 
tio equals 0.015 which is much less than 5.8, placing the 
system of Fig[l2]in the SF regime. We must also empha- 
size that the above ratio is difficult to obtain using VPI, 
as this method does not optimize a parameterized trial 
function, but is rather independent of a trial function. 
Next to this, the VMC method failed to "automatically" 
optimize the trial wave function at this high value of Vq 
in order to reach a ground-state. This is the reason why 
we resorted to a "manual" optimization. 

The upper left inset of Fig. [T2l shows a magnified view 
of the VPI density-overlap between the lattice wells. In 
Fig. [Til on e can see that the tunneling of the BEC be- 
tween the lattice wells is substantially reduced, however, 
upon a careful inspection of Fig.[T2"land the inset, one can 
see that some overlap remains between the lattice wells. 
In fact, Eq.® reveals a global SFF of (23.93 ± 0.06) % 
for that system. That is, superfluidity is still present at 
this high Vq. 



2. Large a c 

Second, we tried to obtain a pure MI by increas- 
ing a c substantially. Therefore, Fig.[l3l displays several 
densities {ri2D{x,y)) for the same system of Fig. [3] with 
k = 1.07T, obtained by increasing a c to extreme values: 
a c = 1.00 [frame (a)]; 1.10 [frame (b)]; 1.50 [frame (c)]; 
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VPI, N = 8, V = 10 




FIG. 13: As in Fig-ITU but for V = 10 and (a) a c = 1.00, (b) 1.10, (c) 1.50, and (d) 2.00, respectively. The a c , k and V are 
in units of a^o, o^o an< ^ f^ha, respectively, where a^o = \/h/(mu}ho)- The x and y are in units of a^o, and {ri2D(x,y)} is in 
units of . 



and 2.00 [frame (d)], respectively. The density-peaks at 
the corners of the CHOCL have the largest amplitudes 
as the particles are repelled away from the center of the 
trap. At this a c , the particles prefer to occupy the cor- 
ners to reduce the repulsive potential energy, but in (a) 
and (b) there is still some small probability for them to 
occupy the wells closer to the center of the trap. It turns 
out that for (a) and (b) there remains BEC overlap (tun- 
neling) between the lattice wells and their global SFF is 
(46.31 ± 0.07) % for a c = 1.00 and (39.82 ± 0.31%) for 
a c = 1.10. The (hid(i)) of Fig-FIBTa) displayed in Fig. 
IT41 reveals a remaining BEC overlap between the wells 
signalling the presence of superfluidity. For even larger 
a c , notice that the central wells are almost vacant in both 
frames (c) and (d) of Fig. [T3l as the extremely strong re- 
pulsion has expelled all the atoms from the central cells. 
However, upon inspecting the {nir>(x)) in Fig. [15] corre- 
sponding to Fig.lTBTd'). one can observe for a c — 2.00 a 
still-existing BEC overlap between the lattice wells via 
the central cell, with a global SFF of (22.46 ± 0.05) %. 

To this end, the main question that remains, then, is 
whether our systems above are only an SF, or a mixture 
of coexisting SF and MI regimes. One can possibly an- 
swer this question by computing the occupancies (Num ) 
of the individual lattice wells [50|, using the counting 
method in Ref. 39]. If (iV( ijfe )) is an integer, then one 
can talk about a MI. In that sense, Table Q] presents 
(N(ijk)) obtained by VPI for the densities in Figs. [TT1 and 



[TBT a) — fd). In general, the (N^j^) are all fractions and 
not integers, and (AT( 000) ), (iV (010) ), and (N (011) ) for 
a c = 1.50 and 2.00, respectively. It might then be pos- 
sible to argue that the empty cells in FigsfTBTc) and (d) 
could constitute a empty (vacuum) MI regime. There- 
fore the CHOCL systems always retain a SF component 
coexistent with MI regimes. 

For a decisive check of the latter possibility, the local 
SFF with respect to the total N, i.e. Eq. (|3"6"|) plus the 
condition (|37p . was computed for the same systems of 
Fig- US] a t a c = 1.50 and 2.00. The reason for taking the 
SFF with respect to N instead of (N^j^)), is because the 
wells (000), (010), and (011) have a very low occupancy 
(N(ijk)) for et c > 1.20; they are close to being empty. 
Therefore, it does not make sense to compute the SFF 
with respect to (Nujf.)) as we did for the lower a c in 
Scc lIIIC2l 



Table HT1 presents results for a c > 1.2 and the same 
lattice wells as in Table H] Viewed from this angle, the 
SFF is negligibly small in all lattice wells. There is a 
nonsignificant SFF of ~ 1% for the well (111) at a c = 
1.50, indicating that some SF has migrated to the corners 
of the CHOCL. The lattice wells with a negligible SFF 
define a local MI phase. In contrast, the global SFF is 
significant. 
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TABLE I: VPI occupancies {Nujk)) of the systems in Figs. Illl and Figs. I13f a-d~). plus an additional system at a c — 1.20 (density 
not shown), at various lattice sites (cells) R„ = (ijk) representative of the whole OL. From left to right: a c is the HC diameter, 
Vo is the OL depth, followed by the cells (000), (010), (011), and (111), respectively. Lengths and energies are in trap units 
a ho — ^h/(mu! ho ) and hoj ho , respectively. 



a c 


V 


(000) 


(010) 


(011) 


(111) 


{a-ho) 


(nu) ho ) 




0.02 


60 


1.382 
±3.169 x 10" 3 


0.536 
±2.449 x 10" 3 


0.222 
±1.797 x 10" 3 


0.058 
±9.103 x 10" 4 


1.00 


10 


0.134 
±1.338 x 10 -3 


0.125 
±1.101 x 10" 3 


0.227 
±1.446 x 10 -3 


0.471 
±1.936 x 10 -3 


1.10 


10 


0.090 
±1.255 x 10" 3 


0.075 
±8.760 x 10" 4 


0.167 
±1.270 x 10" 3 


0.597 
±1.885 x 10" 3 


1.20 


10 


0.015 
±3.723 x 10" 4 


0.016 
±2.732 x 10" 4 


0.064 
±5.350 x 10 -4 


0.787 
±9.054 x 10" 4 


1.50 


10 


6.429 x 10" 4 
±7.495 x 10" 5 


2.063 x 10" 3 
±8.418 x 10" 5 


0.027 
±4.349 x 10" 4 


0.837 
±9.144 x 10" 4 


2.00 


10 


1.386 x 10" 6 
±6.864 x 10 -7 


7.444 x 10~ 5 
±2.485 x 10~ 5 


5.472 x 10~ 3 
±1.585 x 10~ 4 


0.583 
±1.471 x 10" 3 



TABLE II: VPI Superfluid fractions {pa/ p) in individual cells [i.e., with respect to the total N using Eq. (|36[) plus the condition 
(|37|) ] for some of the same systems in Table[I] Lengths and energies are in trap units a; 1Q = \Jh/ (mujj and hujho, respectively. 
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FIG. 14: Integrated VPI ID optical density (n 1D (x)) [Eq.(gIJ] 
of Fig. 113( a) along the a;— axis. The a c , k and Vo are in units 
of aho, a^ and Uojho, respectively, where aho = \Jh/ (muho)- 



FIG. 15: Integrated VPI ID optical density (n 1D (x)) [Eq.flHJ] 
of Fig. I13f d) along the a;— axis. The a c , k and Vo are in units 
of aho, aZn ant ^ fi^ho, respectively, where aho = yfh/ {tnojho). 
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In conclusion, the HC Bose gases in the CHOCL traps 
presented here develop coexisting SF-MI phases at ex- 
tremely large HS radii a c . It also turns out, that one 
must carefully distinguish between global and local su- 
perfluidity when it comes to defining local MI domains. 
The presence of a global SF does not imply the absence 
of a local MI phase. The systems presented here display 
a significant global SFF up to a c = 2.00; but the same 
cannot be stated about the local SFF in the lattice wells. 



F. Momentum density 

In this section, we conclude with a computation of the 
k— space momentum distributions of ground-state densi- 
ties. The goal is to check for the presence of momentum 
states higher than k = 0. According to a discussion by 
Mullin 51], Chester [52|,[53[ proved that there can be no 
BEC in k > states unless there is condensation into the 
k = momentum state. 

The momentum density is calculated by a numerical 
Fourier transform (FT) of the spatial density according 
to 



pFT(k x ,k y ) 

\ r+°° 
4^ 



dy{n 2 D(x,y)) exp(-ik • r), 



(42) 



where r = xi + yj, and k = k x i + fcjj, where i and 
j are unit vectors. Fig. 1161 reveals pFT{k x , k y ) from the 
top frame to the bottom frame, respectively for Figs. [TTJ 
and [T3l a) and (d). One observes that there is always a 
broad central BEC peak surrounded by Bragg peaks sig- 
nalling the presence of k > states in the system. This is 
even the case for the system of FigfT^Td). where no cen- 
tral spatial density peak is present. The pft results are 
another manifestation of the Chester theorem mentioned 
above. 

Next to this, had there been only one broad zero- 
momentum peak, this would have indicated the presence 
of a "pure" MI state [HSl]. The particles are therefore 
not locked in their positions as in a MI and rather display 
a mobility, arising from the hopping from one lattice well 
to another. Further, as a result of trapping, our central 
BEC density is not a sharp function of the momentum 
k. This is in line with the finding of Hen and Rigol 
15], who reported that the k = density peak reveals a 
smoother dependence on k in a trapped system than a 
homogeneous system. 

Our momentum distributions in Fig. [TH] display sim- 
ilar features to Fig. 1(a) of Spielman et al. [46]. The 
diffractive structure is indicative of the presence of a 
SF state (their Fig. 17, top frame). As the system pro- 
gresses into the MI regime, the diffractive structure weak- 
ens (their Fig. 17, middle and bottom frames), as the in- 
tensity of the peaks declines. Accordingly, our results 
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FIG. 16: Momentum densities obtained from Eg. (142)1 . Top 
frame: pp-j^fe^, k y ) for Fig. 1111 middle frame: for Fig. I13f a~); 
bottom frame: for Fig.ED^d). The p-p^ikx, k y ) is in units 
of a\ a , and k x and k y are in units of a^, where a,ho = 



^/h/(mujho)- The a c 
huji lo , respectively. 



k and Vo are in units of a,ho, a h „ and 



of Fig. [TH] are in line with those of Spielman et al. [461 ] 
who indicated that (citing them:) "the diffractive struc- 
ture persists deep into the Mott regime" . However, this 
diffractive structure leads to other manifestations as out- 
lined in the next section. In addition, Gerbier et al. 
[5^ experimentally explored phase coherence of ultracold 
Bose gases trapped in OLs. These authors found that 
phase coherence persisted in the MI phase by studying 
the interference pattern of the density distribution of an 
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expanding 87 Rb BEC released from an OL. The persis- 
tence of this interference in the MI phase was attributed 
to short-range coherence fundamentally attributed to the 
presence of partice-hole pairs. 

According to Yi et al. [ll|, the Bragg peaks in the 
momentum distribution should narrow with a rise of the 
repulsive interactions in the system. However, the mo- 
mentum density in Fig. [16] does not reveal this feature 
since even at extreme repulsion a c ~ O(l) the peaks 
have significant width. This is because the bosons are 
still able to hop between lattice sites to another generat- 
ing a quasimomentum distribution [TTI j . The interactions 
have a finite range a c which significantly reduces the lo- 
calization effect as a c becomes larger. It should be noted, 
that at very large a c , such as a c — 2.00, the only way the 
bosonic HSs can move around is by "turning around" 
each other. 

Yi et al. also noted that an external harmonic trap 
broadens the momentum distribution of the bosons in 
an OL, particularly the condensate part. In fact, this is 
what is observed in Fig. [To! 

G. Connecting with the work of Brouzos et al. |l8l) 

In much relevance to our work, Brouzos et al. ex- 
plored earlier ID few bosonic systems in multiwell traps. 
They chiefly studied the effect of repulsive interactions 
on the distribution of particles in these traps. It was 
found, that for a ID homogeneous multiwell trap with 
hard-wall boundaries and a commensurate filling factor 
(their) v = 1, the wells become equally occupied as the 
repulsion between the bosons rises [their Fig. 2(a)]. The 
particles thus tend to reduce their interaction energy by 
this rearrangement. In the weakly repulsive regime, and 
for a filling factor v < 1, i.e. an incommensurate filling, 
the particles are driven away from the central wells to- 
wards the outer wells and do not equally occupy the wells 
[their Fig. 10(a)]. In analogy to the latter result, we find 
that for an (although) inhomogeneous 3D OL with fill- 
ing N/Nl < 1 (i.e. v = N/Nl), the particles are also 
driven away from the central well towards the outer lat- 
tice wells, particularly the corners of the CHOCL. The 
inset in their Fig. 10(b) shows how the populations of 
their wells change with the repulsion between the bosons. 
Most importantly, whereas the central site (s = 4) gets 
less populated as their g rises, the site (s = 1) farther 
away gets more populated. Our Figs. [3] and [5] display a 
similar behavior for the 3D CHOCL systems, although 
there does not seem to be a transfer of particles in the 
intermediate lattice wells between the central well of the 
CHOCL and its wells at the corners, as they reported. 
This is because none of the intermediate wells gains and 
subsequently loses population as those demonstrated by 
the inset of their Fig. 10(b). 

In contrast to their ID inhomogeneous result with 
v > 1 [their Fig. 14(a)], where even at large g the central 
wells remains significantly populated, the central cell of 



the inhomogeneous 3D CHOCL becomes almost vacant! 
Thus, there is a higher tendency for the central well to 
become vacant in a 3D CHOCL than a ID multiwell trap 
with harmonic confinement. As such, the dimensionality 
of this system seems to play a crucial role in this spe- 
cific feature. This is particularly since the particles in 
3D do not face the "energy obstacles" between the wells 
as much as in a ID multiwell trap with harmonic con- 
finement p^| . when the tunnel from the center to the 
corners of the CHOCL. These "energetic obstacles" force 
the particles in a ID inhomogeneous OL with v > 1 to 
favor the occupation of the middle wells. 

For the filling factor v < 1, their particles remained 
delocalized and coherence did not vanish completely, con- 
trary to the commensurate filling case. Our CHOCL sys- 
tems with v < 1 are also incommensurately filled and as 
such will always display a derealization of particles and 
an associated coherence manifested by a remaining SFF. 

With regards to their momentum distributions in Fig. 
11(c), they display a structure with Bragg peaks, where 
the central peak is lowered because of a partial loss of 
coherence with increasing repulsion between the bosons. 
Similarly, our Fig. 1161 also displays a lowering of the cen- 
tral peak with increasing repulsion. The Bragg peaks 
surrounding the central one in our Fig. [TBI reveal accord- 
ingly incomplete localization and persisting coherence, 
even as the system enters into the MI regime. In the 
work of Brouzos et at, the Bragg peaks of the momen- 
tum distribution of a uniform ID Bose gas in a multiwell 
trap with v — 1 vanish into a broad, smooth, Gaussian 
peak centered at k — with a rise of g [their Fig. 5(a)]. 

Most importantly, Brouzos et al. computed the CF in 
the lowest natural orbital (£ — 0) [their Fig. 11(b)] and 
showed how it declines with increasing g, whereas it rises 
for higher ones (£ > 0). Our Figs. [3] and [5] display the 
same result for I = 0, where for (000), (001), and (010) 
the CF decreases. However, it rises for (011) and (111). 

Our results for the 3D CHOCL with v ^ 1 are then 
very much in line with those of Brouzos et al. [18| for the 
ID inhomogeneous case with v ^ 1. 



IV. CONCLUSIONS 

In summary, then, we have presented a numerical in- 
vestigation of the CF and SFF of a few HC bosons in a 
3D CHOCL of 3 x 3 x 3 lattice sites. The global and local 
CF and SFF were computed for each CHOCL well and 
were explored as functions of a c and the lattice wave vec- 
tor k — ir/d, d being the lattice spacing. The role of the 
interference between the condensates in all lattice wells, 
and its effect on the CF in one well was also studied. In a 
major part of this work, an achievement of a "pure" MI 
state was attempted, yet only mixed SF-MI phases were 
obtained. 

The most important result of this paper is an oppos- 
ing behavior for the global CF and SFF as functions of k. 
Whereas the global CF increases with increasing k, the 
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global SFF decreases with k. For the CF, this was ex- 
plained by approximating each lattice well by a HO trap 
and the preference of the bosons to occupy the lowest HO 
level as the energy-level spacing increases with increasing 
k. For the SF, this was explained by an increase in the 
localization of the bosons in each lattice well as the local 
confinement strength increases with increasing k. 

From a local perspective, it was found that the con- 
densate is depleted with a rise of a c in the lattice wells 
R„ = (000), (001), and (010), whereas it rises with in- 
creasing a c for (011) and (111)- This is because of the 
tunneling of the condensate away from the trap center to 
the farther lattice sites at the edges of the trap. The local 
CF is also enhanced with a rise of k in all lattice wells due 
to the same reasons outlined for the global CF. Further, 
when interference with all-neighbor sites is included, the 
CF in each lattice well is enhanced beyond the case with 
no interference effects. The local SFF as a function of a c 
in each lattice well displays a similar behavior to the CF 
as it overall declines in the wells (000) and (010) whereas 
it rises in (011) and (111). In contrast to the local CF, 
however, the local SFF decreases with increasing k, re- 
vealing again an opposing behavior to the CF on a local 
scale as well. An important point to mention, is that the 
local CFs and SFFs are of the same order of magnitude. 

The energy of the system rises with an increase in a c , 
largely because the average onsite interactions rise at a 
rate overwhelming the drop in kinetic energy, or mobil- 
ity of the bosons. Further, the energy does not tend to 
stabilize towards a fermionization limit as it happens in 
the ID homogeneous multiwell trap with commensurate 
filling i/ = l@. 

Finally, the possibility for achieving a "pure" MI state 



in a 3D CHOCL with filling v < 1 was investigated in 
order to compare with the ID homogeneous multiwell 
trap with hard-wall boundaries explored by Brouzos et 
al. [l8j]. This was performed by either increasing Vq to 
60 (huiho) or by increasing a c to huge values a c ~ O(l). 
It was found that no "pure" MI could be achieved in a 3D 
CHOCL, but rather a coexisting SF-(vacuum)MI phase. 
This is because for a filling v < 1, the particles remain 
delocalized and coherence cannot vanish completely even 
at huge a c . From another point of view, the small number 
of 3 x 3 x 3 sites does not furnish the ground for a flat MI 
density [l3T ] to appear in the system. Consequently, our 
systems remain superfluid, even in the extreme repulsive 
regime with a global SFF of - 20% for a c > 2. The 
presence of a global SFF does not guarantee the absence 
of local MI regimes, however. Further, the empty lattice 
wells in the CHOCL define vacuum MI regimes. 

The momentum density of our systems with v < 1 
reveals a broad central BEC peak surrounded by Bragg 
peaks, even at a c ~ 0(1). This proves again the presence 
of a coexisting SF-MI phase in the extremely repulsive 
regime due to the derealization of the particles. Thus, 
our results are in line with those presented by Brouzos et 
al. [lH . Had there been only one central BEC peak (2|| , 
the systems would have been defined as a "pure" MI. 
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